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Abstract 

We describe a systematic approach to generating approximate actions for 
the lattice simulation of QCD. Three different tuning conditions are defined 
to match approximate with true actions, and it is shown that these three 
conditions become equivalent when the approximate and true actions are suf- 
ficiently close. We present a detailed study of approximate actions in the 
lattice Schwinger model together with an exploratory study of full QCD at un- 
physical parameter values. We find that the technicalities of the approximate 
action approach work quite well. However, very delicate tuning is necessary 
to find an approximate action which gives good predictions for all physical 
observables. Our best view of the immediate applicability of the methods we 
describe is to allow high statistics studies of particular physical observables 
after a low statistics full fermion simulation has been used to prepare the 
stage. 

12.38.Gc, 11.15.Ha, 02.70.Lq 
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I. INTRODUCTION 



The design of improved algorithms for simulation with dynamical fermions has provided 
an important challenge for lattice QCD. Perturbative arguments and numerical experiments 
suggest that a significant component of the effect of dynamical quarks can be accounted for 
by a shift in the effective lattice spacing. The quenched approximation takes this view to its 
extreme, and ignores all dynamical quark effects except for this shift in scale. In most cases, 
the systematic errors introduced by quenching have so far turned out to be of the same order 
as the best presently achievable statistical errors in full dynamical simulations. As available 
computer power grows, and as we study larger volumes and smaller lattice spacings, we 
can expect this situation to change, and to begin to see significant discrepancies between 
quenched and full theory results. Systematic errors due to quenching will then turn out 
be larger than full theory statistical errors, and the quenched approximation will loose its 
usefulness. 

In a recent paper it has been shown that the quenched algorithm is a single mem- 
ber of a large class of approximate QCD simulation algorithms. The members of this class 
interpolate in a smooth fashion between the extreme cases of quenched QCD simulations 
on the one hand, and full dynamical QCD simulations on the other. Individual members 
can be considered as different approximations to the full dynamical QCD action. So long 
as the quenched approximation represents the only significant example of an approximate 
simulation algorithm, it is possible to consider it in isolation and to accept its historically 
important contributions to lattice QCD while remaining wary of its mathematical correct- 
ness. However, once one begins to consider variants which seek to improve on the quenched 
approximation without themselves being exact, a number of important issues arise, and a 
more critical approach is needed. 

First, of course, is the question of what we mean mathematically by the term 'approxi- 
mate action'. The qualitative view of quenched QCD which we normally adopt is that, at 
fixed lattice spacing and coupling, it produces ensembles which are reasonably close to the 
ensembles which would be generated for full QCD at some shifted coupling. Reasonable 
closeness is however a vague term. What it tends to mean in practice is that ensemble 
expectations of a set of operators in quenched and full QCD seem to agree within statistical 
errors, so long as one allows for an effective lattice spacing shift. 

In this paper we consider the general principles which can be applied to define what is 
meant by a simulation algorithm based on an approximate action in the sense defined by the 
quenched approximation. We will describe a systematic approach to generating approximate 
actions and will discuss conditions which allow the tuning of any candidate approximate 
action to optimize its matching to a true action of interest. The methods we develop, while 
motivated by the approximate action view, are also applicable whenever we wish to match 
results from different lattice simulations. In particular, the tuning we propose defines the 
optimum matching conditions between Wilson and Kogut-Susskind fermion simulations, and 
between naive and improved action simulations. 

A second issue of concern, when considering simulations with approximate actions, is the 
question of their computational performance. The results reported in were obtained on 
very small systems, and for these systems it was found that approximate action simulations 
were very competitive with exact simulations including full fermions. However, a question 
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left unanswered was whether this competitiveness would continue when the method was ap- 
plied to more realistic situations. The basic proposal of [Q was to attempt an approximation 
of the fermion determinant in QCD using Wilson loop operators of increasing size. On a 
small system there are relatively few such operators, but the number of operators grows 
very fast as the size of the system increases and the computational work required would also 
grow. Some preliminary data which we reported in [0 for the case of the Schwinger model 
seems to bear out this view. There we found that only very slow improvements in the quality 
of the approximation were achieved as Wilson loops of systematically increasing size were 
added to the approximate action. An independent study of QCD at realistic lattice sizes 
0] was also somewhat negative. This study found that the plaquette operator was actually 
a very poor approximation to the fermion determinant in typical large simulation. Both of 
these results suggest that very many Wilson loops would be needed to get a numerically 
useful approximation in realistic situations. 

We report in this paper, however, that the pessimistic view just described is an artifact of 
the particular approximation prescription defined in The systematic approach to tuning 
approximate actions which we develop here allows us to generate much more aggressive 
approximations than are otherwise possible, and we find that it is not necessary to include 
Wilson loops of all shapes and sizes to get useful approximate actions. Indeed, the best 
approximate actions in the Schwinger model and QCD cases which we have studied are 
found to include only a few Wilson loop operators of quite surprising shape and dimension. 

An issue which we do not address in this paper is that of the behaviour of approximate 
actions in the continuum limit. The quenched case is clearly problematic in this limit [Q. 
However, the strategy we propose here is a numerical one in which one may choose the 
most appropriate approximate action for any given simulation at fixed lattice spacing and 
coupling. If one varies the lattice spacing and coupling, then the corresponding approximate 
action would also vary. It seems conceivable to us that the variations in approximate action 
which result would be sufficient to allow a mathematically correct continuum limit of the 
approximate action to be taken. 

The plan of the remainder of this paper is as follows. In Section |T|, we define a simple 
cumulant expansion method which allows direct comparison of the expectations of operators 
calculated in two different path integral measures. We argue that this expansion method 
is the most efficient way of determining the shifts introduced in expectations of physical 
observables due to changes in the underlying path integral measure. 

The concept of an approximate action is defined in Section |T|. The quenched QCD action 
is the prototypical example here. Approximate actions allow us to define corresponding 
approximate simulation algorithms, and we develop in this section three different ways in 
which approximate actions might be defined. We then argue that all three definitions become 
equivalent when the approximate action is sufficiently close to the action which it attempts 
to approximate. 

As an initial test, we consider the Schwinger model. Section |^ gives the explicit form 
of this two-dimensional model which we have used, and describes some necessary technical 
features of the calculations. Simulation results are presented in the following two sections. 
In Section [V|, we consider the quenched approximation in some detail. Physical quantities 
calculated here include the static potential and the mass of the mesonic bound state of the 
model. We present comparisons of these observables calculated in the quenched approxi- 
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mate simulation and full dynamical simulation. Section ^ describes a variety of different 
extensions to the quenched approximation. The analysis begins with the scenario proposed 
in in which the fermion determinant is expanded in terms of Wilson loops of increasing 
size. We show first that this expansion technique is quite straightforward to execute and 
appears to be numerically stable. We argue, however, that this ordered expansion is very 
slow to converge. Instead, we propose an aggressive approach to developing approximate 
actions which directly targets the interesting physics. This approach is seen to improve 
considerably the slow convergence apparent in 

In Section |V11| , we describe a trial application of the approximate action approach to 
QCD. The purpose of the test is to gain some understanding of the potential of the method 
for QCD simulation and to identify practical problems which remain to be solved. The test 
described is for a small lattice with heavy quarks. 

Finally, our conclusions are presented in Section [VIII . 



II. THE COMPARISON OF PATH INTEGRAL MEASURES 

In all lattice simulation, the goal is to calculate expectation values of interesting operators 
in some path integral measure. The measure is given by some action which we think of as 
the true action and denote by St- Our goal is to find an action Sa which allows us to 
approximate the measure generated by St. The prototypical example is given by lattice 
QCD. The fundamental degrees of freedom here are the SU{3) link matrices U^{x) assigned 
to edges which join nearest neighbouring sites on a four dimensional hypercube. A possible 
'true' action for this system is the lattice QCD action which includes a Wilson fermion 
determinant term, 

St = -I3Y1 ^ReTrf/n - Tr log M+M, (2.1) 

while a possible approximate action is the quenched lattice QCD action, 

'^a = -/3'E^ReTrf/n. (2.2) 
□ 'J 

Un here represents the product of gauge links around a fundamental plaquette, while is 
the Wilson fermion matrix for hopping parameter k. We have also implicitly allowed for a 
bare coupling shift between exact and approximate cases by including different /9's in the 
different actions. 

Expectation values of operators F in the path integral measures generated by St and Sa 
are denoted respectively by {F)t and (-F)a, and are given by 

Before addressing the problem of how to choose an approximate action, we first need a 
technique to compare these expectations. The simplest numerical comparison technique is 
to generate independent sequences of configurations for each measure, to evaluate expec- 
tations on these sequences, and then compare. For our purposes, this simple approach is 
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not computationally feasible. We are attempting to find approximate actions for a given 
true action because the computational work needed to generate independent configurations 
in the true measure is very large. Our eventual goal is to find a tuning technique which 
allows us to generate an approximate action without ever having to simulate the measure 
of the true action. Thus we are lead to consider methods by which we can use ensembles of 
configurations generated in one measure to estimate expectations of operators in a second 
measure. 

Provided the fundamental degrees of freedom of two action functionals, 5*1 and 5*2, are 
the same {i.e. provided the configuration space upon which they are defined is the same), 
we can directly relate expectations in the corresponding measures, 

= nSr^ = + «^ - ^^^^^^^^^ ' ^2.4) 

where 

A,, = S,- S, = -A,, (2.5) 

and 

QiA) = j^e\ (2.6) 

This comparison formula is exact. However, the calculation of the overlap expectations of 
Ci{Aij) with F — {F)i is very difficult numerically, since Aij is typically a quantity which is 
of 0{V), where V is the volume of our system. If we use this formula, therefore, to calculate 
expectations for Sj by correcting an ensemble of configurations generated for Si, we find 
that we must re- weight with factors which are O(e^). The fluctuations in expectations so 
calculated are also of O(e^), so the work needed to get reasonable statistics is exponential in 
V. As a result, a direct application of this formula to compare {F)i and {F)j is not feasible 
numerically. 

We adopt instead a cumulant expansion technique which is generated by expanding 
Cj(A) in a Taylor series, 

Q{A) = 1 + (A - (A),) 

+1((A-(A),)^ + (A),2-(A2),)+---. (2.7) 

Substituting this expansion in (|2.4|), we find 

(F),-(F), = ((F-(F),) (A,, - (A,,),)), 

+ l{iF-{Fh) (A,, - (A,,),)^), + ■ ■ ■ . 

(2.8) 

The leading term on the right hand side here is the connected part of the correlation of the 
operator F with the difference in the actions Aj^. As before, Aij is a quantity of 0{V)] 
but now we don't have to exponentiate, and so we expect fluctuations only of 0{V) rather 
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than O(e^). Since any numerical simulation must be carried out at finite volume, we can 
expect that at least this leading correction between the expectation values in the different 
measures will be calculable with finite work. Of course, higher order terms will require 
calculation of connected correlations involving higher powers of Aj^. These correlations will 
have correspondingly larger fluctuations. The connected correlations involving A^- will have 
fluctuations of 0(y"), for example, and will be correspondingly harder to calculate. One 
should view this cumulant expansion technique, therefore, as the numerical version of an 
asymptotic expansion. If the difference A^- is sufficiently small, one or perhaps two terms 
in this expansion may be calculable with reasonable work. 



III. TUNING THE APPROXIMATE ACTION 

In our prototypical example introduced in the last section, the true action St is the 
lattice QCD action which combines a single plaquette gauge term and a Wilson fermion 
determinant term expressed as the trace of the log of the Wilson hopping matrix. Our 
prototypical approximate action Sa includes just the single plaquette gauge term but with 
a possibly different coupling from that in the true action. More generally, we imagine that 
the approximate action depends on a finite number of couplings 

Sa = Sa{ai, . . . ,an) . (3.1) 

One possibility, for example, is that Sa includes a sum of different translationally invariant 
Wilson loop operators [|T|. The are then the couplings multiplying each different Wilson 
loop operator in the sum. Another possibility is to replace the fermion determinant term with 
an approximate bosonic term proposed by Liischer 0. The Oj then include the plaquette 
coupling for the gauge term together with all the roots of the polynomial used to approximate 
the inverse of the fermion matrix in the Liischer approach. The question we now address is 
how to choose the couplings so that Sa becomes a good approximation to the true action 
St- We define three different tuning schemes for the coefficients, and show that these three 
different schemes are closely related. 



A. Minimizing the Distance Norm 

The first tuning procedure is described in There it is observed that the measure 
defined by an action functional S allows us to define an inner product ( , ) on the space of 
gauge invariant functionals of the link matrices, which takes the form, 

(F, G) -{FG)- —^j^UY^z^ , (3-2) 

where F, G, and S are all functions of the link matrices U^{x). This inner product turns the 
space of gauge invariant functions into a Hilbert space, and allows us to define the distance 
d{F, G) between any two functions as 

d{F, G) = {{F - G), (F - G)) = {{F - - G)) . (3.3) 
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Functions of link matrices can now be considered as vectors in this Hilbert space and the 
natural tuning condition in this language is to minimize the distance, or equivalently to 
maximize the overlap, between the two vectors which represent the true and approximate 
actions. Since an overall constant shift in an action functional has no effect on the measure it 
generates, the tuning should minimize the connected part only. Thus, the first prescription 
is to minimize 

d{Sa - {Sa),St - {St)) = {{Aat - {Aat) f) (3.4) 

over the space of couplings ai, . . . ,a„. Aat is the difference in the approximate and true 
actions, 

Aat = Sa-St = -Ata (3.5) 

and the object to be minimized is simply the variance a'^{Aat) of this difference, 

a\A,t) = ((A„i - (A,,))2) (3.6) 

evaluated in the measure ( ) . Differentiating with respect to each of these couplings in turn, 
we generate a set of n simultaneous equations to determine the minimum, 

{{S, - {S,)){Aat - (Aat))) = for z = 1, . . . , n, (3.7) 

where 

d 

Si = — — Sa{ai, . . . , a„) . (3.8) 

In deriving this particular tuning condition, we have been purposefully vague about the mea- 
sure in which expectations are to be evaluated. There are two obvious measures available, 
that generated by 5*^, and that generated by Sa- Either measure will allow a calculation of 
the tuned values for the couplings Oj, so this first tuning condition actually has two possible 
variants. Our preferred choice is to minimize in the true measure. This minimization is 
an explicit problem, but in principle requires generating a set of configurations in the full 
measure. In practice, we could also use the cumulative expansion methods of the last section 
to estimate expectations in the true measure using our best guess as to the approximate 
measure. Minimizing in the approximate measure, on the other hand, is an implicit prob- 
lem since we only learn the correct values for the coefficients after the minimization is 
completed. Minimization in this case requires an iterative method which first guesses the 
correct couplings, then evaluates the corrections to those guesses. 

For notational simplicity in what follows, we will assume that equations of the form ( p.7[ ) 
apply for all values of i, and we will drop the phrase 'for i = 1, . . . , n' in such equations. 



B. Operator Matching 

A different approach to tuning the approximate action is suggested by the observation 
that, in some cases at least, both true and approximate actions represent different regu- 
larizations of the same continuum theory. For the prototypical case of full Wilson fermion 
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action versus quenched action this is, of course, not correct since the quenched approximate 
action does not have a unitary continuum hmit. However, it is possible to imagine cases 
where this is indeed true. For example, one could consider the case where one action was 
the full Wilson fermion action, while the second action was an 0(a) improved action. In 
this view, the couplings in both actions are bare parameters in different regularizations of 
the same continuum theory. 

Normally, one fixes the mapping from bare to renormalized couplings by specifying renor- 
malization conditions for each bare parameter whose finite part needs to be fixed. These 
conditions directly relate regularized theories to the renormalized continuum theories which 
they regulate. In our case, imposing renormalization conditions which relate the regular- 
ized theories to continuum renormalized theories is not feasible, since this requires that we 
separately generate a full continuum limit extrapolation of each regularized theory. What 
one can do, however, is to impose renormalization conditions directly between the two lat- 
tice regularized theories. If these conditions are maintained in the passage to the continuum 
limit, then both regularized theories will proceed to the same renormalized continuum theory 
as we require. 

Consider, therefore, matching the approximate action Sa{ Cli , . . . , (In ) with the true action 
Sf. The approximate action has n couplings which must be determined, so we must impose 
n independent matching conditions between the approximate and true measures to fix these 
n couplings. The obvious possibilities include matching operators such as m^, nip, rriN 0. 
These particular operators however are not ideal for our purposes since they are difficult 
to measure accurately on a small number of configurations, and they are not particularly 
sensitive to small changes in the couplings which we need to determine. The choice 
of matching conditions which we make, however, is not really very important since any 
set of conditions will correctly fix the cutoff dependent divergences of the couplings being 
matched, and will differ only in the way they fix finite parts. Different matchings will differ 
by finite renormalizations only, so we are free to choose some simply implemented matching 
conditions. An alternative set which is simple to calculate, and which should be sensitive to 
changes in the Oj, is to match the operators Si = dS/dai between the approximate and true 
actions. 

Thus, the second tuning condition we propose is to match the expectations of the n 
operators Si {i = 1, ... ,n) in the approximate and true measures, 

{S,)t = {Si)a. (3.9) 

Applying the cumulant expansion to {Si) a, we find 

{S,)t = {S,)t + {{S, - {S,)t)Cti-Aat))t , (3.10) 

which directly gives 

{{Si-{S,)t)Ct{-Aat))t = 0. (3.11) 

C. Maximizing Acceptance 

Consider using an approximate action as part of an exact Markov process which generates 
configurations distributed according to the measure of the true action. The particular exact 
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algorithm we propose involves a standard Markov transition step which satisfies detailed 
balance for the approximate action. This step generates new trial configurations which do 
not have the correct true action distribution. One can then correct to the true distribution by 
executing a Metropolis accept/reject on the trial configurations. Our third tuning condition 
for the approximate action is to maximize acceptance in this exact algorithm. 

Let Pa{U U') be a Markov transition probability for link configuration U to go 
(reversibly) to link configuration U'. If this transition probability satisfies detailed balance 
for the approximate action Sa, then 



-Sa 



U 



PaiU ^U') = e-""^ PaiU' ^U) 



(3.12) 



where the notation indicates that the relevant action terms are evaluated on the link 
configurations labeled by U and U' . To get an exact algorithm for the true action St we must 
now execute an accept/reject step. The probability, Aat, to accept a trial step generated by 
Pa must satisfy a modified detailed balance, 



Aat{U ^ U') 



,Aat{U' ^U) 



and the optimum choice for this acceptance probability is 

Aat{U U') = min (l, e^u^vi^at) 



where 



(A 



at) 



A 



at\u' 



A 



at\u 



(3.13) 



(3.14) 



(3.15) 



Maximizing acceptance for this compound algorithm requires maximizing the weighted av- 
erage of Aat over all configurations U which can occur at the start of a trial step, and all 
configurations U' which can occur at the end of such a step. If we denote this weighted 
average as {{Aat)), we have 



{{Aat)) = / d[U']d[U] A 



at 



Pa{U ^ U') 



U 



Jd[U]e-s^ 

Since Pa satisfies detailed balance for action Sa, it is straightforward to show that 
which implies 

{{Su^U'iAat))) < 0. 



(3.16) 



(3.17) 



(3.18) 



Further, if {{6{Aat))) is small, and if only the leading terms in the Taylor expansion of 



are significant, the acceptance probability is given approximately by [0 



{{Aat)) 



erfc - 
\2 



-((^(AaO)) 



(3.19) 



Acceptance, therefore, is maximized when —{{6{Aat))) is minimized. The worst case for 
acceptance occurs when the transition probability Pa{U U') defines a perfect heat bath 
transition for the approximate action. 
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PaiU ^ U') 



e ^° 



Jd[U]e- 



(3.20) 



Trial configurations U' in this case are completely independent of the configurations U which 
they attempt to replace. In this worst case, we have: 

-((5(A.,))) = -((Aa*|^,- Aatl^)) 
= -{{Aat)a- {Kt)t) 

= -{{Aat - {Aat)t)Cti-Aat))t ■ (3.21) 

Maximizing acceptance requires maximizing this last term over the couplings oi, . . . , in 
the approximate action Sa- If we differentiate with respect to the couplings Oi, we find the 
following condition: 

((5, - (5,)i)a(-A„t))i - ((5, - {S,Cti-A,t))t)iAat - {A,t)t)Cti-Aat))t = 0. (3.22) 

An important question to ask once an approximate action has been tuned to maximize 
acceptance is whether the acceptance achieved by this maximization is large enough to 
generate a practical exact algorithm. The answer, of course, is given by the value of {{Aat)) 
which results. {{Aat)) is given in ( |3.19| ). We can rearrange this equation by using the worst 



case formula for —{{6{Aat))) given in (|3.21| ). When expanded, this latter equation takes the 
very simple form, 

-{mat))) = {{Aat-{Aat)ty)t+o{Ai) 

= a^{Aat) + 0{Al,). (3.23) 

(The subscript t on denotes that the expectations involved are to be evaluated in the 
true action measure). Thus, we find 

{{Aat)) ^ erfc ^y-{{6{Aat)))) ^ erfc (^a,(A„,)) • (3.24) 

The actual minimum value for {{Aat)) which defines a practical algorithm is a function of 
the ratio of work involved in generating uncorrelated trial configurations with approximate 
and true actions. In QCD, the generation of uncorrelated configurations with, for example, 
an exact hybrid algorithm is very expensive. This suggests that even very small accep- 
tances might suffice to generate practical exact algorithms based on approximate actions. 
A 50% acceptance is achieved when (A^t) ~ 1, while a 10% acceptance is achieved when 
a^{Aat) ^ 5. 



D. Tuning Procedure 



We have, at this point, proposed three different tuning conditions, (|3.7|, 3.11, p.22|), to 



define the adjustable parameters in an approximate action. Two of these equations involve 
the cumulant expansion factor Ct(±Aaj). As we have argued in the last section, this factor 
is not directly calculable in a single numerical experiment with a single ensemble generated 
by either the approximate or true action. Instead, to apply these tuning conditions, we must 
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expand Ct{±Aat) in powers of Aat using ( |2.7|) . The tuning conditions which result when we 
keep just the leading and next to leading terms of this expansion take the following forms. 
First, if vector overlap in the true measure is maximized using ( p.7| ), we have unchanged 



■at 



{A,t)t))t = 0. 



(3.25) 



Second, if operators are matched in approximate and true measures then we have 



{{S, - {S,)t){Aat - {A,t)t))t - l^S, - {S,)t){Aat - {Aat)tY)t + ■ ■ ■ = . (3.26) 



Finally, if acceptance is maximized using ( ^.22] ), we have 



((5, - {S,)t){Aat - {Aat)t))t - 7((5. - {S^)t){Aat - {Aat)t f)t + " " " = . (3.27) 



These three conditions all agree to leading order in Aat, and only begin to disagree at next 
to leading order in the cumulant expansion. If we have managed to find approximate and 
true actions which agree sufficiently closely that the cumulant expansion is a numerically 
useful tool for comparing expectations in the two measures, then we would expect the terms 
at next to leading order in these equations to be much smaller than the leading order terms. 

These observations suggests the following self-consistent approach to generating approx- 
imate actions for lattice simulation. First, apply one of the three conditions listed here 
to tune the effective action. This tuning can be carried out for an ensemble of configura- 
tions generated either for the approximate or for the full measure, so long as the cumulant 
expansion is used appropriately. Once the approximate action is determined, a variety of 
operators should be calculated in this approximate action, and further, the shift between 
expectations in approximate and true actions should be determined for these operators us- 
ing the cumulant expansion. If the shift in operators so calculated is small, then we expect 
that the approximate action provides a good representation of the true action. Further, we 
expect that all three possible tuning conditions will be approximately satisfied, and that 
the cumulant expansion will work well for most operators. On the other hand, if the shift 
in the operators so calculated is large, then the approximate action is not a good represen- 
tation of the true action. In this case, we find that the different matching conditions will 
produce significantly different values for the approximate action couplings. The situation 
which will then prevail is that we will be able to adjust certain operators to match between 
the approximate and true measures. However, this adjustment will not cause other different 
operators to match, and simulation with the approximate action will not provide a useful 
approximate ensemble for the study of properties of the true action ensemble. 

At leading order in the cumulant expansion, the common tuning condition from all three 
methods takes the form 



As written, this condition requires that expectations are to be evaluated in the true mea- 
sure. As we suggested above, we can use the cumulant expansion to apply this condition 
using operators evaluated in the approximate measure also, so long as we use the cumulant 
expansion to correct appropriately. We have 



((5. -(5.),) (A, 



■at 



{Aat)t))t = 0. 



(3.28) 
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(3.29) 



These two forms differ only by terms which we have aheady neglected when we work to 
leading order. Thus we can quite consistently choose to drop the next to leading terms on 
the right hand side of ( |3.29| ) when we wish to use an approximate action ensemble rather 
than a true action ensemble to perform tuning. 

Our leading order approximate action tuning condition now takes the simple form 

{iS,-{Si)){Aat-{Kt))) = (3.30) 

where expectations are evaluated in either true or approximate measures. This condition 
represents a set of coupled equations for the tunable coefficients in the approximate 
action Sa, and are precisely the equations which minimize a'^{Aat) over these coefficients. 
cr^(Aaf) therefore provides a key measure of the quality of an approximate action. It also, 
according to ( p.23|) , determines whether a practical exact algorithm can be built using a 
global accept/reject step to correct from approximate to true action. In what follows, we 
use this variance to discriminate between candidate approximate actions and to guide our 
search for an exact update algorithm based on the approximate action. 



E. Wilson Loop Approximations 

To close this section, let us consider the technicalities of minimizing cr^. The approximate 
actions which we study in the remainder of this paper are simple sums of Wilson loop 
operators of different sizes and shapes. We will use the notation Wj to denote a particular 
Wilson loop operator which we always take to be translationally and rotationally invariant. 
A typical example of such an operator in QCD is the 1x1 plaquette summed over all sites 
and orientations, 

Wixi = E^ReTrf/n. (3.31) 

□ 3 

Wilson loop operators are always summed over all distinct orientations and locations on 
the lattice. Thus the operators we consider are rotationally and translationally invariant. 
This is sufficient for our purposes since our major effort is to find an approximation for the 
fermion determinant which is also rotationally and translationally invariant. 

The class of approximate actions which we will consider now takes the simple form 

n 

S'a(ai, . . . ,a„) = -^ajWi, (3.32) 

i=l 

and our goal will be to find the best approximation within this class to the true action 

St = -(3W,^,-T (3.33) 

where T, defined as. 
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T =^Tr log M^M (3.34) 

is the trace of the log of the fermion determinant for Uf equal mass fermions. We have 
immediately that 

d 

—Sa{ai,...,an)=Si = -Wi. (3.35) 
The difference between approximate and true actions for this case is 

^at = Sa — St 

i 

= T-Y,a[W, (3.36) 
1 

where, for convenience, we have absorbed the pure gauge term by redefining the coefficients 

a'lxi = ctixi - /3 (3.37) 
a[ = ai i^lxl. (3.38) 

and the work of finding an approximate action reduces to that of finding an approximation 
for T of the form T ^ J2i cay^i- For convenience in what follows, we drop the prime super- 
script on the parameters a, which define this approximation. Our lowest order minimization 
condition ( |3.30| ) is now simply expressed as 

n 

- m)m - (>v,)))«, = m - mm - (m . (3.39) 

To determine the parameters Oj which define the best approximation to T we must calculate 
connected correlations of each of the Wilson loop operators included in the approximate 
action both with themselves {i.e. ((WH>Vi))(W~(Wj)))), and with T («.e. {iWi-{Wi)){T- 
(T)))). Once these correlators are available, then the evaluation of the a^'s is achieved by 
solving the linear system ( |3.39| ). 

It is also useful, in the process of generating the coefficients of the approximation to T, 
to define a Gram-Schmidt orthogonalized basis of Wilson loop operators. If the original loops 
are Wj for i = 1, . . . ,n, then the orthogonalized loops are VVj and are defined recursively by 

= Wi- < Wi > (3.40) 

and 

i-l 

W, = W~ <W,> t>l, (3.41) 

where 



13 



C-». (3.42) 

The operators Wj therefore satisfy 

(>V,W,) = {W^)Si, , (W,) = . (3.43) 



and, in terms of these orthogonahzed operators, the solution of the hnear system ( p.39| ) is 
given by 

a.^<^-(^-<^>)>. (3.44) 



IV. THE SCHWINGER MODEL 

In order to test the ideas developed in the proceeding section, we need a model which 
exhibits the important features of the QCD lattice model, but is computationally accessible. 
The Schwinger model offers a suitable testing ground. This theory is a two-dimensional gauge 
theory with gauge group U{1) rather than SU{3). With no fermions present the theory is 
(almost) trivial. All the interesting dynamics are generated by fermions, and we expect very 
significant dynamical fermion effects to be present as a result. Computationally, it is possible 
to generate full dynamical fermion simulations in just a few days on a workstation, so the 
model allows considerable experimentation with multiple different actions over reasonable 
time scales. The fact that much is known exactly and perturbatively about the corresponding 
continuum theory provides an important bonus. 

The explicit form of Schwinger Model action which we adopt is 

S=-(3Y.^eUn - ^Tr log MtM,. (4.1) 

Uf is the number of equal mass fermions included, which we take to be 2 in all our simulations. 
Gauge elements f/^(x) are defined on the links of a two dimensional square lattice and take 
values in the group U{1). The hopping matrix is given by the usual Wilson form 

M,(x, y) = 6^^y - kK{x, y), (4.2) 
K{x, y) = ((1 - j,)U,{x)6,+^,y + (1 + j,)Ul{x - A)4-a..) • (4-3) 

M can be decomposed into its red/black (even/odd) components. 



M^1-kK= ' 7 . (4.4) 




Using this decomposition, one may construct a Hermitian, positive definite, red/black pre- 
conditioned matrix H^^^ such that 

H^^) = M(")M(^)t , M(^) = 1 - K^K^"K"^ . (4.5) 
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For two flavours of equal mass fermions, the fermion contribution to the Schwinger model 
action is then given by 

T = Tr log MlM^ = Tr log . (4.6) 

In what is to follow, we will use preconditioned iterative methods to calculate inverses and 
logs of the fermion matrix. We have found that the red-black form here provides the best 
preconditioning . 

In the language of the last section, our true action is the Schwinger model action with 
dynamical fermions included and is given by ( |4.1| ). We used the hybrid molecular dynamics 
algorithm with a flve step integration scheme [§ to generate ensembles of conflgurations 
distributed according to this true action. All the approximate actions which we considered 
involved varying sets of Wilson loop operators. To generate ensembles distributed according 
to these actions, we used a simple Metropolis link update algorithm. 

The only other simulation issue worth commenting upon is that of the calculation of 
T and for individual conflgurations in an ensemble. These calculations are required 
whenever we need to flnd correlations between action terms. T is deflned in ( |4.6| ) as the 
trace of the log of the matrix H^'^K To calculate T we adopt a stochastic approach A 
number A^^ of Gaussian distributed vectors (pai^) are generated. The index x here runs 
over the (even) sites of the lattice, and a = 1, . . . , A^^ labels the different vectors generated. 
We use the notation { )^ to denote expectations over the vectors and the probability 
distribution of the components of these vectors are normalized so 

Mx)(l)liy))^ = 6a,b5x,y ■ (4.7) 

The stochastic estimators, and Erp2, which we adopted for T and respectively, are 
given by 

Et = ^Y: UiL<Pa) (4.8) 



a 



where, for convenience, we deflne 



i-- Y: (0l^0a) (j>lLcl>,) (4.9) 



L = logH^^K (4.10) 



We have immediately that 



{Et)^ = TtL = T (4.11) 
(Et2)<^ = (TrL)2 = (4.12) 

((ET)^),= (TrL)2 + -|-Tr(L2). (4.13) 

Note that Et^ as deflned here is an unbiased estimator of T^. The other obvious choice of 
estimator for is {EtY. We see from these equations that E"^ is a biased estimator for 
which only becomes unbiased in the limit A"^ oo. The variances of the estimators Et 
and Erp2 are also easy to calculate. We have 
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(4.14) 



The variance in Et^ has a more comphcated form, but a similar l/N^ overall factor. When 
using these estimators to calculate T and T^, we must therefore choose iV<^ sufficiently large 
that the fluctuations due to are small. 

The remaining technical problem is to calculate Lcp = log H^^^cf). We use a Chebychev 
polynomial approximation (of order A'",:;) to perform this calculation. To optimize the poly- 
nomial approximation used, it is necessary to know the maximum and minimum eigenvalues 
of the matrix H^^\ These values are obtained with a Lanczos method The calculation of 
the Chebychev polynomial approximation itself then requires Nc matrix-times-vector prod- 
ucts for each vector 0. In all our analyses, we have found that the time needed to estimate 
T and is a significant fraction (5%-50%) of the time needed to perform a full dynamical 
simulation. 

Arbitrary numbers of equal mass fermions can, of course, be dealt with by incorporating 
a factor nf/2 along with T in the accompanying formalism. 



V. LOWEST ORDER RESULTS 
A. The Quenched Approximation 



Our preliminary application of the methods introduced in Section |T| to the Schwinger 
model was performed on lattices of size 16^ and 32^, at P values from 2.0 to 3.0, and k 
from 0.25 to 0.265. Typically, we generated 1000 or 2000 configurations using either the 
true full fermion action or an approximate action. Maximum and minimum eigenvalues of 
(e) y^QYQ calculated on each of these configurations using a Lanczos algorithm. The major 
computational effort is then to evaluate the stochastic estimators Et and E^^ of T and 
respectively on each configuration. We used 30-50 different Gaussian 0's per configuration 
to accumulate these estimates, and we used Chebychev polynomials of order 100-200 to 
calculate L(f) = log H^^^(j) for each 0. We normally worked with configurations which were 
first fixed to Landau gauge. This was expected to reduce the infiuence of gauge artifacts in 
the Lanczos procedure and allowed the use of Fourier acceleration methods. In practice, we 
found that preconditioning additional to that provided by the red/black construction gave 
little extra benefit. 

We also accumulated measurements of a variety of different Wilson loop operators on 



each configuration. As described in Section [III D| , we will use the notation Wj for such a 



loop operator labeled by the generic index i. Where appropriate, the generic index will 
be replaced by a descriptive index. For example, Wixi denotes the simple 1x1 plaquette 
Wilson loop operator which, for the Schwinger model, is defined by 

Wixi = El^ef/n. (5.1) 
□ 

Typical loop operators considered include 1x1, 2x1, 2x2, 3x1 etc. 

When the raw data for the estimators Et, Ej-2 and the loop operators Wj have been 
calculated on a configuration by configuration basis, the final analysis required is to evaluate 
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correlations of different loop operators with Et, and to derive the values of the approximate 
action coefficients according to the conditions of Section [IIID . 



The first and simplest case to consider is the quenched approximation. The true action 
(for two flavours of equal mass fermions, ( [4.1|) ) is 

St = -(3J2 Ref/n - Tr log M^M = -/5Wixi - T . (5.2) 

□ 

The approximate action which defines the quenched approximation includes just the single 
plaquette Wilson loop operator, Wixi- The parameters of this approximate action must 
however be adjusted according to the conditions in Section [111 \J\ and so we write it with a 
different coupling, /?'. Thus, 

Sa = -/3'Wixi (5.3) 
The difference between these two actions is 

A = Aat = Sa-St = T + <5/3Wixi (5.4) 

where 

Sf3 = f]-f3'. (5.5) 

6P is the shift in plaquette coupling induced by approximating T with a plaquette term only. 
Our lowest order tuning condition is that the variance, (T^(A), is to be minimized in either 
the true or approximate gauge measure. Explicitly we have 

a\A) = (((T - (T)) + 5/3(Wixi - (Wixi)))') • (5.6) 

The condition on 5(3 to minimize this expectation is 

((Wixi-(Wixi))(T-(T))) 
^" ((Wixi-(Wixi))^) ■ 

In the initial test, we used the quenched algorithm to generate configurations. Thus the 
ensemble which we have generated is the appropriate one for the approximate action, and we 
choose to minimize the variance of A in this measure. For (3' = 2.50 and n = 0.26, we found 
5/3 = — 0.212±0.010. These and all other errors were evaluated using the bootstrap method, 
with subensembling to detect and remove autocorrelations. The interpretation of this result 
is that, within the space of quenched actions defined by (|5.3| ), the best approximation to 
the true action ( ^.2|) with /? = 2.500 — 0.212(±0.010) and k = .26 is given by the quenched 
action with p' = 2.500. 

Having defined an approximate action, it is immediately important to check the quality 
of the approximation which that action generates. In particular, we ask whether we have 
been justified in using just the leading term in the cumulant expansion when tuning the 
approximate action. To investigate this, we generated a new ensemble of 2000 configurations 
using the hybrid molecular dynamics algorithm to incorporate dynamical fermions {k = 0.26) 
a.t [3 = (3' + 6/3 = 2.29. As a simple preliminary test, we compared (Wixi) as measured 
using the approximate and true actions. The quenched, approximate action value was 
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(Wixi)a = 783.3 ±0.3 13' = 2.50 (5.8) 

while the true dynamical fermion action value was 

(W^^i)^ = 782.8 ±0.3 ±1.4 /3 = 2.29, k = 0.26 . (5.9) 

The error ±1.4 in this last equation is due to the uncertainty with which we have deter- 
mined 6(3. The agreement between the two measurements is a very encouraging result. The 
expectations of Wixi match in the true and approximate measures. We are therefore in the 
situation where tuning using the common lowest order tuning conditions ( 3.30| ) produces an 



approximate action which matches the expectation of an operator in the approximate and 
true measures. This is our first evidence that, in this instance, we have been justified in 
taking only the leading term in the cumulant expansion. 



B. Bound state mass 

We now investigate how well the quenched approximation performs in calculating ex- 
pectations of some interesting physical operators in the Schwinger model. The continuum 
model with a single species of massless fermions has a 'vector' bound state with mass e/ 
where e is the fermion charge. For n/ flavours of massless 'quarks', there is an SU{nf) 
symmetry, leading to an iso-singlet 'vector' particle with mass 



m^ = eJ^ (5.10) 
V vr 

and a massless {nj — l)-plet [|l^. On the lattice we have two flavours and look for a 



bound state by studying the large Euclidean time dependence of zero momentum correlators 
involving tpjitp. The latter were constructed from 'quark' propagators which were obtained 
by a Conjugate Gradient solver applied to the red/black preconditioned Hermitian matrix 
^(e) QY by the BiCGStab algorithm |jll[] used with M^^\ The latter turned out to be the 



more efficient method. 

Figure |I] shows the vector state effective mass m„ obtained from a full dynamical fermion 
simulation with the true action at (3 = 2.29, k = 0.26 compared with that obtained from a 
quenched simulation at = 2.50. The quenched values plotted for rriy include the leading 
'deficit correction' which is defined by the first term on the right hand side of the cumulant 
expansion formula ( |2.^ ). The dashed line in this figure serves both to indicate the (uncor- 
rected) value for obtained by a quenched simulation with (3' = 2.29, and the Euclidean 
time range used to extract mass estimates. From the true action simulation and from the 
deficit corrected estimate from the quenched simulation respectively, we find = 0.369(3) 
and rriy = 0.367(7) where the bracketed digit is the error in the last significant digit. These 
vales are in complete agreement. In the deficit corrected value, the contribution due to the 
j3' shift from j3' = 2.29 to (3' = 2.5 is —.041 and that due to the deficit correction term is 
-.037. 

This is pleasing but not yet surprising. The approximate action simulation is again 
generating predictions which are in agreement with those from the correct true action simu- 
lation. It must, however, be pointed out that the statistical fiuctuations in the approximate 
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action estimates are, for the same number of configurations, somewhat larger (x2.5). So 
any comparison of algorithm efficiency must take this into account. We will return to this 
issue later. 

The relatively small difference between quenched and full theory estimates in this partic- 
ular observable is in accord with experience gained in lattice studies of the hadron spectrum. 
It also means that there is little point in using m„ when studying more complicated approx- 
imate actions. To really stretch the method we now study an observable which is explicitly 
dependent on the long distance dynamics. 



C. Static potential 

Figure ^ shows a comparison of the static potential using quenched simulation data at 
(3' = 2.50, full dynamical simulation data at /3 = 2.29, k = 0.26 and the deficit corrected 
estimate from the quenched data. The results have been expressed in dimensionless units us- 
ing our result m„ = 0.37 in lattice units. Also shown is the (exact) infinite volume quenched 
result (dashed line) and the continuum massless Schwinger result (full line) obtained by 
evaluating Wilson loops due to a vector field of mass m^. The continuum result is 

^^^ = ^(l-e— «). (5.11) 

Note that, expressed this way, the continuum result is independent of the number of fiavours. 
The comparison between the lattice estimates is satisfactory although the errors associated 
with the deficit estimate ( |2.8[ ) get quite large at large distances. This could signal a system- 
atic under-estimate of fermion effects at increasing R, but it is not yet statistically significant. 
The same behaviour is seen in a direct comparison using Wilson loops with increasing area 
(see below). It seems therefore, that the potential screening effects due to light dynamical 
fermions are capable of being reproduced by these first order estimates based on quenched 
configurations. 



VI. EXTENSIONS 

A. Naive Improvement 

The results of the last section are quite encouraging. Approximating T with a 1 x 1 
plaquette is seen to shift significantly the quenched results towards the full theory results. 
The absolute measure of how well we have actually done was defined in the discussion 
following ( p.23| ) and is given by the corresponding values of ct^(A). For the 32^ lattice 



at p = 2.50, K = 0.26, we find ^^(r) = 99.0(2), while a^{T + 6pWi^i) = 3.1(7). These 



values allow us to predict acceptances according to (|3.24|) for algorithms which first generate 
trial configurations, then accept/reject those configurations according to the change in T or 
T + 6pWixi respectively. The value of o"^(T) = 99.0(2) implies that an algorithm which uses 
the pure gauge term of the Schwinger model at unshifted (3 to generate trial configurations 
will have an acceptance of 2 x lO'^^ (erf c( 0.5^991))). The value (j2(T + 5/3Wixi) = 3.1(7), 
on the other hand, implies that an algorithm which uses the pure gauge term at a properly 
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tuned shifted (3 will have an acceptance of 0.2. For the Schwinger model, at the parameters 
we have analyzed, the 1x1 loop is therefore seen to be a reasonable approximation to the 
fermion determinant, since it produces a finite (non zero) value of acceptance. It is not 
perfect, however, since acceptance is still quite small. 

We wish to consider now whether, and to what extent, the addition of further loops 
to the approximate action can improve the situation. For this and all remaining analysis 
in this paper, we will concentrate on the effects of approximate actions on loop operators 
only. The single loop analysis just described shows that the loop expectation values are 
the most sensitive to changes in action. The effective mass is rather insensitive, and in any 
case, is already quite well described by the addition of just the 1x1 plaquette term in the 
approximate action. 

The first addition which we considered involves an approximate action where we included 
loops containing up to three plaquettes. The loops involved are the 1 x 1, 1 x 2, and 1x3 
rectangles (denoted by 1, 2, and 3 for short), together with a 'chair' shaped loop of length 
8 (denoted by 3'). These four loops are shown in Figure |^ and are some, but not all, of the 
simple loops of perimeter length 8 or less. An immediate problem with this loop expansion 
is that the loop possibilities expand very rapidly as the loop length increases. Even at this 
simple expansion level, we are already ignoring some non-trivial length 8 loops (for example, 
the 1x1 squared loop operator). 

The target full fermion action of our analysis has (3 = 2.50 and k = 0.26. Expressed in 
terms of Gram-Schmidt orthogonalized loop operators, the generic approximation of T is 
now 

T = aiWi + asWs + asWs + as'Ws' • (6.1) 

Initially we have no information about the values of the coefficients a^, so we set them to 
zero and simulate with an approximate action which contains only the pure gauge part of 
the Schwinger action with /3 = 2.50 {i.e. we simulate with the unshifted quenched version 
of the Schwinger action). The ensemble so generated allows us to make new estimates of 
the coefficients a^, which can be used to generate a new approximation. We repeated this 
process four times in all. Table | summarizes the sequence followed. 

Ensembles of 2000 configurations were generated for each approximation, and the ap- 
proximate action coefficients for the next approximation were calculated on these ensembles 
using the formulae of Section [ill E| . The results achieved are given in Table 0. Approxima- 



tion in a quenched simulation with (3 = 2.50. The value which results for af*^ from this 
approximation is 0.212(10) which is exactly the value of —6(3 determined in the last section. 
Here however, we have specified that the true action has (3 = 2.50. Thus we interpret af''' 
as defining the best quenched approximation coupling shift to use to simulate the full the- 
ory with j3 = 2.50. This best quenched theory has (3' = 2.50 + 0.212(10) and defines our 
second approximation (approximation 1 in Table ||). Each approximation of course gener- 
ates a different ensemble, and the values for the approximate action coefficients Oj therefore 
change from approximation to approximation. This is clear also in the table. But a very 
nice feature to note is the fundamental stability of the procedure we have adopted. The co- 
efficients which result after approximation 3 are, within errors, the same as the coefficients 
after approximation 0. 
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For the highest approximation attempted (approximation 3), we also hst the correspond- 
ing root variance o"(Wj) and the product aj(j(Wj). We have, when the coefficients Oj are 
tuned according to p.44| , that 

a\T - <^^y^^) = ^\T) - Y: a,V(>V.) . (6.2) 

i i 

Thus, ajCr(VVi) is a measure of the importance of a particular loop operator in approximating 
T. If the addition of loops of increasing size generated a rapidily converging approximation 
to the true action, we would expect ajO"(>Vj) to decrease rapidly as the loop size increases. 
It is clear from the table, however, that this is not the case for the Schwinger model. 
The contributions of two and three plaquette loops are not significantly smaller than the 
contribution of the single plaquette loop. 

In Figure ^, we show the static potential, at fixed distance in lattice units (for R = a 
and R = 5a), versus the approximation number, n. We show both the raw measurements 
made using each set of configurations, and the deficit-corrected estimates given by ( |2.8| ) as 
described above. The corresponding dynamical fermion result (2,000 configurations) is also 
shown for each inter-quark distance. As expected, we see that adding loops of increasing 
size has very little effect on the short distance potential aX R = a. However, at i? = 5a 
there is an effect, and a decreasing trend is evident in the (uncorrected) data obtained as 
larger loops are added to approximate action. This decrease is quite slow unfortunately, 
and con&ms that the systematic addition of loops of increasing size produces only a slowly 
converging approximation to the true action. 

One significant positive observation can made from Figure ^ This observation, antic- 
ipated in the previous section, is that at each approximation order, n > 0, the deficit- 
corrected estimate given by ( p.8|) is compatible with the dynamical fermion result. At ffist 
sight, it may therefore seem that there is no need to go beyond n = 1, i.e. to use other than 
quenched data to make properly corrected estimates of the full theory. However, as empha- 



sized in Section one can only be confident of such cumulant expansion-based estimates 
when the shifts involved are known to be small. This will only be so when the corrections 
are based on a suitably accurate approximate action. 

Note finally that the corrections are not simply shifts in /?, but are non-trivial effects 
of dynamical fermions. One may think of the estimates as a form of re-weighting in which 
the approximate vacuum used becomes closer to the true vacuum as more loops are added. 
However, what is very clear is that the addition of small loops only in the approximate action 
improves physics only at the scales appropriate for small loops. Thus the static potential 
at short range is almost perfectly determined by the approximate action. At longer range, 
however, the approximate action is not very successful in improving the quenched result. 



B. Selective Improvement 

This analysis suggests that an approximate action based on a systematic expansion of T 
in terms of loops of increasing size is going to require the inclusion of a very large number 
before convergence is achieved. However, the impression already gained is that, once a 
loop of a particular size is added to the approximate action, the addition of more loops 
of a similar size is not particularly productive. This suggests an alternative approach of 
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adding loops selectively rather than systematically. The aim would be to include typical 
loops of different sizes rather than all loops up to a given size. As an experiment, we chose 
to consider approximate actions with 1 x 1, 3 x 3 and 5x5 loops included in turn. For 
simplicity, tuned coefficients defining all these actions were determined on a single quenched, 
(3 = 2.50 ensemble. 

The results achieved by these new approximate actions are given in Table |T|. The left- 
most column shows the values for various loop expectations from a full dynamical fermion 
simulation. The remaining columns gives the values achieved by the various different approx- 
imate actions. What we see is that the expectation value achieved for a loop in a simulation 
which includes that loop in the approximate action correctly matches its full djTiamical sim- 
ulation value. This is in good agreement with our analysis of Section |TTI| where we argued 
that our approximate action tuning conditions are such as to match loops included in the 
approximate action to their correct values in the true action. Loops not explicitly included 
in the approximate action are not so well tracked. There is however considerable improve- 
ment as we add more loops to the approximate action. For example, the 3x8 loop changes 
from value 3.7 when only the 1x1 term is included in the approximate action to the value 
10.5 when 1x1,3x3, and 5x5 loops are included. This is to be compared with its exact 
full theory value of 14.4. 

A check on the results achieved is to execute an iteration of the tuning algorithm to see 
if our results are sensitive to the choice of coefficients in the 1x1 + 3x3 + 5x5 approximate 
action used. To execute this iteration, we used our Ix 1 + 3x3 + 5x5 approximate action 
to generate a new ensemble on which we recalculated the coefficients 01x1,03x3 and 05x5. 
Table |IV| shows the results of this retuning. Again we find that the results are quite stable. 
In the Schwinger model, a set of quenched configurations seems sufficient to calculate the 
approximate action coefficients. 

C. Optimized Improvement of The Effective Action 

As a final test of the method, we attempt a very direct approach to defining the effective 
action by using an ensemble of configurations generated using the full fermion algorithm. 
The ensemble is for a 16^ lattice at /5 = 2.50 and k = 0.26. In all we generated 1,000 config- 
urations, and systematically considered a large number of different possible approximations 
for T of the form 

T^^a^W, (6.3) 

i 

where the index i here runs over some selected subset of the Wilson loop operators 
which we have calculated. For example, one possible approximation to consider is T ^ 
O2X2VV2X2 + Oix3Wix3- For each possible approximation, we choose values of Oj which min- 
imize (j^(T — J2i c^iVVi) for that approximation, and then determine the approximate action 
which generates the minimum value for cx^. The results of this minimization are listed in 
Table 0, and the optimum loop choice is shown in Figure ^. The first line of the table gives 
the value of cr^(T) with no approximation attempted. The second line indicates that the 
best single loop approximation to T is obtained when the loop is the 3x3 loop. The best 
two loop approximation is obtained by using the 3x3 and 1x6 loops, etc. The measure 
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of the quality of a given approximation is the value of cr^(T — J^iCa^i) listed in the final 
column of this table. The best one loop approximation reduces from 2.69 to 1.02. Adding 
more loops reduces cr^ further, but not at the same rate as adding only a single loop. 



D. Stability 

In all cases considered so far, we have tuned the effective action by minimizing cr^. In Sec- 



tion m we proposed three different tuning conditions. These all agree to lowest order in the 
cumulant expansion. They differ at higher order, however. To test whether the differences 
between the three tuning conditions are significant, we compared the coefficients for the 



effective actions tuned according to (|3.71 ), and (|3.22|) . The results of this comparison 



are shown in table |VI[ The values determined by minimizing for the coefficients of the 
3x3, 1x6, 5x5 and 1x2 loops in the best four loop approximation for T are given in the 
second column of the table. If instead, we tune to match the expectations of the four loops 
in the approximate action between the true and approximate measures, then the values of 
the coefficients are shifted by the amounts listed in the third column of the table. If we tune 
to maximize acceptance, then the shifts necessary are given in the final column of the table. 
In all cases the shifts are very small, compatible with zero at this level of statistics. This is 
indicative that we have found an approximate action which is sufficiently close to the true 
action that all the tuning conditions which we have considered are self-consistent. 

Having arrived at a tuned approximate action [e.g. coefficients given in Table |V|), we 
have made preliminary tests of an exact algorithm as described in Section [III C\ A similar 
exact algorithm was introduced in the context of the bosonic algorithm by Peardon 



and further developed by others [|13[. These tests have shown that correct results can indeed 
be obtained with an experimental acceptance as predicted in Section |III C|. However, careful 



tuning is required when applying the algorithm to large systems. The approximate action 
must yield a^{Aat not too different from 1 so that reasonable acceptances are obtained. 



VII. APPLICATION TO QCD 

So far, the main numerical results in this paper have concerned the application of ap- 
proximate actions to the Schwinger model. We chose this model because it allowed useful 
comparisons between the ensembles generated by an exact full fermion algorithm and those 
generated by a variety of different approximate actions in reasonable computational time. 
However, our main interest in the application of approximate actions is in lattice QCD. A 
complete repetition of the analysis so far described for lattice QCD is a considerable compu- 
tational undertaking, and beyond the scope of this exploratory paper. Nevertheless, it has 
been possible to gain some indication of both the potential of the method and the practical 
problems which still remain for lattice QCD. 

The exploratory simulation which we have carried out is for full QCD on a small lattice. 
The true action which we considered was the standard single plaquette pure gauge action 
with coupling (3, combined with a two-flavour Wilson fermion action with hopping parameter 
K. The lattice size adopted was 6^ with f3 = 5.7 and k = 0.16. This is obviously a very 
small lattice with very heavy quarks. The basic goal was to determine the extent to which 
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it is possible to approximate the true action for this theory with an action made up from a 
finite set of Wilson loops. 

The direct approach to analyzing this question follows exactly the strategy adopted at 
the end of the last section, and is simply to generate an ensemble of configurations using 
the full QCD action, to evaluate T and various different loops on these configurations, and 
then to determine the set of loops which best approximates T. This is the approach we have 
adopted. We used hybrid molecular dynamics to generate configurations. The ensemble 
analyzed contained 516 configurations. Each configuration was separated by 10 trajectories 
of length 0.5. T was calculated as for the Schwinger model on each configuration. We 
also calculated a variety of different Wilson loop operators on each configuration. These 
loops were chosen to be representative rather than complete. Four dimensions allow many 
more possible loop configurations than two dimensions, and an exhaustive study of all the 
loop possibilities is not justified in the present instance where we are working with a very 
unphysical set of lattice parameters. Our goal is simply to explore the possibilities rather 
than to present a definitive calculation. 

Following our experience with the Schwinger model, we looked at a variety of different 
loop shapes and sizes. These included two, three, and four dimensional loops, with total 
length varying from 4 to 12 lattice spacings. The only non-trivial loop of length 4, of course, 
is the pure gauge 1x1 plaquette action, while the length 12 loops considered included 
5 X 1, 4 X 2, and 3x3 2-dimensional loops together with more complicated three and 
four dimensional constructions. A particular loop construction is defined by specifying the 
steps in the various possible directions which a link in that loop can represent. Thus the 
specification {1212} denotes a loop generated by first stepping in the direction 1 with positive 
sense, then in direction 2 with positive sense, then in direction 1 with negative sense, and 
finally in direction 2 with negative sense. This set of four steps brings the loop operator 
back to its starting site. We can then take the trace in the usual way. The digits 1 and 
2 can represent any two distinct directions in the four dimensional lattice. An unbarred 
digit defines a step in the positive direction along the corresponding axis, a barred digit 
defines a step in the negative direction. The specification implicitly assumes a sum over all 
independent orientations and starting locations of the loop. For example, the specification 
we have just written represents the simple 1x1 plaquette summed over all orientations and 
lattice sites, and is therefore the standard one plaquette pure gauge action term. 

To test how well the true action can be approximated, we then seek to minimize a'^{T — 
J2i diWi) over the different loops considered. In all, we calculated over 40 loops. This allows 
us to consider a very large number of different approximate actions. These possibilitys 
include actions with a single loop include, actions with two loops included etc. 

Table |VII| shows some of the results we achieved. The first column lists the variety of 
loops which we considered. The order in which these loops are listed represents their relative 
importance in reducing fluctuations in the quantity cr^(A) and is determined only after all 
analysis is complete. The second column lists the minimum achieved for cr^(A) when just a 
single loop is included in the approximate action. This minimum, of course, varies depending 
on which particular loop is considered as the single term approximate action, and the second 
column shows this variation. The raw value of the variance in T is listed at the end of the 
table for comparison purposes. 

It is immediately noteworthy that there is a large variation in the values achieved for 
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cr^(A), indicating that there is a large variation in the effectiveness with which different 
loops can individually approximate T. The effective loops tend to be loops with four dimen- 
sional structure, and it is quite striking that the simplest loop of all, the single plaquette 
loop (defined as {1212}) is actually quite unimportant. The best match between true and 
approximate actions is given by the loop which is the first entry in the table. This loop, 
on its own, reduces the fluctuations in St from a value of 295(20) to a value 8.7(4.4). This 
is a remarkable 97% reduction in the raw fluctuations of St. The loop which does this is 
not a simple loop but rather an oddly shaped four dimensional construction which has total 
length 12. We have not been able to identify any remarkable features about this particular 
loop relative to all the other loops which we considered. Other similarly shaped loops are 
seen to reduce the fluctuations in St almost as much, so it is reasonable to imagine that the 
important feature is the four dimensional nature of the loop and its approximate size rather 
than the exact details of how the loop is generated. 

The determination of the best single loop approximate action is only the first analysis 
step. The next is to determine the best two loop approximate action. Our approach is 
exactly as before. We minimize cr^(A) for all two loop approximate actions. In fact, we only 
need to minimize o"^(A) over loops with 1 and n where n is any of the other available loops. 
This minimization determines that the approximate action which uses the loops placed 1 
and 2 in the table is the best two loop approximate action. The final column of Table |V11| 
indicates the value which results for cr^(A) from this minimization. The minimization for 
higher numbers of terms in the effective action can now proceed by iteration. The final 
result we find is that the best n loop approximate action is the action which contains loops 
1, 2, . . . , n. The values which result for o"^(A) with this best n term approximate action are 
all given by the final column of Table |VI1] . 

The second interesting feature of this exploratory study is that the addition of extra loops 
beyond the first loop in the approximate action is not particularly helpful. The approximate 
action with a single loop reduces the raw fluctuations in St from 295(20) to 8.7(4.4). Adding a 
second and third term reduces these fluctuations only to 7.2(4.1) and 5.95(4.00) respectively. 
These shifts are very much smaller than that achieved by a single term alone. 

As even more loops are added, the error in o"^(A) becomes greater than its actual value, 
and we conclude that we have insufficient data to calculate reliably the best approximate 
action for more than a small number of terms, so the actual ordering after the first three 
or four terms should be considered suspect. However, a significant observation which can 
be made is that the 1x1 plaquette loop contribution comes at position 35, which is very 
far down the list. We conclude therefore that the pure gauge plaquette action is actually a 
rather poor approximant of the the full fermion action. 

VIII. CONCLUSIONS 

We have proposed a general framework for the construction and evaluation of approx- 
imate actions for lattice QCD. We have demonstrated how the coefficients in the simplest 
form of such actions can be evaluated directly from a true action simulation or by using 
an iterative procedure which begins with an approximate action simulation, and never re- 
quires simulation with the true action. We have also shown that the results of the iterative 
procedure are quite stable, and converge very quickly. 
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As a first application, we have tested approximate actions based on the systematic ex- 
pansion of the trace of the log of the fermion matrix in Wilson loops of increasing size as 
originally proposed in . This expansion works well for observables sensitive to short range 
physics. Using the Schwinger model, we were able to test longer range observables where it 
becomes clear that an impractically large number of loops will be required to study physics 
over a range of scales using small lattices spacings. 

The systematic error correction estimate proposed by |I| is shown to give reliable esti- 
mates of observables based on approximate action simulations. The better the approximate 
action, the more reliable these 'deficit' estimates become. However, since these estimated 
corrections require a good, unbiased, stochastic estimate of T, the trace of the log of the 
fermion matrix, their practical value is dependent on further progress in developing efficient 
calculational methods for this stochastic estimation. The computational methods which we 
have used here to evaluate T are quite demanding. Typically they require work equivalent to 
that needed to evaluate all meson and hadron propagators at a number of different hopping 
parameters on a given configuration. 

Having established the above behaviour, we have gone on to show that the greatest deter- 
minant of the quality of an effective action is not the number of loops included, but the choice 
of scales represented by the loops which are included. We have found dramatic evidence of 
this both for the two-dimensional test model and for QCD itself. The systematic matching 
and tuning methods have allowed us to propose approximate actions significantly closer to 
the true QCD action than the quenched approximation or, indeed, of an approximation 
based on the ordered inclusion of loops of increasing size. 

We presented two methods for capitalising on approximate actions in lattice simulation. 
The first, as originally proposed in |[^, is to simulate with the approximate action and to 
make (presumably) reliable estimates of the small corrections to desired observables. The 
challenge to progress in this strategy is to improve stochastic evaluation methods for T. The 
second method described is to use the approximate action to generate candidate update 
configurations within an exact algorithm, and we proved by construction that such an exact 
simulation was feasible in the Schwinger model example. The challenge to the construction 
of such an exact algorithm in QCD is to implement the fine tuning required to achieve viable 
acceptance rates. It is perhaps too optimistic to imagine that this fine tuning can be achieved 
by considering actions made of Wilson loop operators only. The QCD test described here 
suggests that the return achieved in adding more Wilson loops to the approximate action 
becomes very poor after the first few loops. On the other hand, the tuning methods we 
have defined are completely general, and can be applied to other candidates for approximate 
actions. In particular, a combination approximate action including a small number of Wilson 
loop operators and a multiboson Liischer ^ term seems to us to be an interesting possibility 
as a tunable candidate approximate action. 

While awaiting these developments, an immediate and exciting application of the meth- 
ods we have described is to perform relatively high statistics measurements using the ap- 
proximate action directly. The strategy here is to make low statistics measurements of the 
quantity of interest using the full theory. These can then be used to help tune an approxi- 
mate action to reproduce this measurement approximately (and other more accessible ones 
to higher levels of accuracy). The approximate action can then be used to give high statistics 
estimates of the quantity of interest {e.g. glueball mass). 
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TABLES 



TABLE L Sequence of approximations to T used to test the effects of adding loops of in- 
creasing size to generate an approximate Schwinger model action. The loops included in these 
approximations are defined in Figure |3|. The superscripts (n) on the coefficients a-"^ indicate the 
approximation used to calculate those coefficients. Note that, in all cases, approximation k uses 
coefficients determined by simulating with approximation k — 1. 



Approximation Form of T used 

T = 

1 T = af^Wi 

2 T = a'i^Wi + a^^^Ws 



TABLE IL Numerical values for the coefficients ai in the expansion of T = Yli O'i^i ^-iid their 
relative contributions (see text). 





(0) 


«1" 


(2) 


(3) 




aP-fOV.) 




0.212(10) 


0.208(10) 


0.197(11) 


0.222(10) 


10.2(2) 


2.26(11) 


m 


0.085(7) 


0.083(7) 


0.073(7) 


0.083(7) 


13.8(3) 


1.15(10) 


m 


0.041(6) 


0.018(11) 


0.034(7) 


0.048(8) 


12.7(2) 


0.61(10) 




0.025(4) 


0.026(5) 


0.034(5) 


0.029(4) 


23.9(4) 


0.69(10) 



TABLE in. Rectangular Wilson loop measurements made on configurations generated by full 
dynamical fermion simulation (true), by an approximate action with 1x1 loop added (column 2), 
with 1x1 and 3x3 loops added (column 3), and with 1 x 1, 3 x 3, and 5x5 loops added (column 
4). In a given column, expectations of loops included in the corresponding approximate action are 
marked with an asterisk. 





true 


1 X 1 


+3x3 


+5x5 


1 X 1 


805.1(3) 


805.5(2)* 


804.9(2)* 


804.4(2)* 


2x2 


415.7(9) 


391.8(7) 


410.4(8) 


408.2(8) 


3x3 


158.9(11) 


118.9(9) 


162.4(10)* 


160.6(10)* 


4x4 


51.2(10) 


23.7(7) 


47.0(8) 


46.7(9) 


5x5 


15.1(8) 


3.8(7) 


8.9(7) 


13.8(6)* 


1 X 2 


640.2(6) 


633.5(4) 


637.6(4) 


636.4(4) 


1 X 3 


512.2(8) 


498.2(5) 


507.5(5) 


506.3(6) 


1 X 4 


411.2(10) 


392.0(6) 


403.9(6) 


402.9(7) 


1 X 5 


330.9(11) 


308.5(7) 


321.5(7) 


320.8(8) 


1 X 8 


172.9(11) 


150.4(7) 


162.0(8) 


162.3(9) 


2x6 


88.8(9) 


58.4(6) 


77.9(7) 


78.0(8) 


3x8 


14.4(6) 


3.7(5) 


10.0(5) 


10.5(6) 
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TABLE IV. Comparison of Wilson loop expectations calculated with a true full dynamical 
simulation, and with approximate actions containing 1 x 1, 3 x 3 and 5x5 loops. Column 2 
lists the loop expectations determined with approximate action coefficients tuned on a quenched 
ensemble. The approximate action coefficients in column 4 are tuned on an ensemble generated 
from the approximate action of column 2, and represent therefore a retuning of the coefficients 
generating column 2. Columns 3 and 5 list the corresponding deficit corrected loop expectations 
for each approximate action. 



loop 


true action 


quenched tuning 


quenched tuning 
(corrected) 


retuned 


retuned 
(corrected) 


1x1 


201.5(02) 


201.2(02) 


201.2(02) 


201.5(01) 


201.5(01) 


2x2 


105.2(06) 


101.8(08) 


104.4(11) 


103.3(05) 


104.3(08) 


3x3 


40.5(07) 


39.8(07) 


39.8(07) 


42.6(08) 


42.6(08) 


4x4 


12.9(06) 


12.4(05) 


13.1(10) 


13.7(03) 


14.1(08) 


5x5 


3.6(04) 


3.7(04) 


3.7(04) 


5.1(04) 


5.1(04) 


1 X 2 


160.7(04) 


159.2(04) 


160.2(06) 


159.8(02) 


160.2(03) 


1 X 3 


128.9(05) 


126.5(06) 


128.5(07) 


127.5(02) 


128.5(04) 


1 X 4 


103.4(06) 


100.5(07) 


103.2(09) 


101.6(03) 


102.9(06) 


1 X 5 


82.9(06) 


79.8(08) 


83.2(12) 


81.0(03) 


82.4(08) 


1 X 8 


42.5(07) 


39.5(08) 


43.2(14) 


41.1(05) 


42.9(12) 


2x6 


22.2(06) 


19.2(07) 


21.0(13) 


20.5(03) 


21.6(06) 


3x8 


3.8(03) 


2.8(04) 


2.0(08) 


3.6(04) 


3.3(08) 


5x7 


1.4(03) 


0.8(02) 


2.0(06) 


0.9(02) 


1.2(06) 


5x8 


0.4(02) 


0.4(02) 


0.4(05) 


-0.1(03) 


1.2(05) 



TABLE V. The residual variance corresponding to an approximate action incorporating the 

Wilson loop terms shown. Th(^ \"ariauc(' is m measure of the (|ualily of the ai")i")roximale act ion. 



Terms Loops a^<(r -Y^iO-i^i) 

2.69 (27) 

1 3 x 3 1.02 (09) 

2 +1 x 6 0.61 (07) 

3 +5 x 5 0.59 (07) 

4 +1 x 2 0.31 (05) 
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TABLE VI. Effect of tuning condition choice on the loop coefficients in a tuned approximate 
action. Loop coefficient values as determined by minimizing arc given in the second column. 



1 ri o trlniTT"Q 
-L lie ollllLo 


in these values when higher order tuning 


coiitiiiioiio aie iiicul poi aieu. 


dlt: y,lVt:ll 111 LlltJ 


remaining 


columns. 








coefficient 






Operators 


Acceptance 






Minimized 


Matched 


Minimized 


a3x3 




0.2150 (12) 


-0.013 (45) 


-0.026 (94) 


^6x1 




0.0325 (18) 


-0.001 (5) 


-0.001 (10) 


05x5 




0.0105 (20) 


-0.004 (5) 


-0.006 (8) 


aix2 




0.0109 (22) 


-0.003 (9) 


-0.004 (18) 
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TABLE VII. Approximate action optimization data for the lattice QCD test described in the 
text. Wilson loop operators are defined in column 2. Column 3 lists the value of cr^(T — ajWi) (no 
sum on i) for the loop operator defined in that row. Column 4 lists the value of a^{T — X)a;=i o,iWi) 
where the sum includes all the loop operators at that row or higher in the table. The order in 
which loop operators are listed is the order which minimizes o^{T — YL\=i when at most i 

loops are included. Note that we present information only on those loops which generate better 
approximate actions than the single plaquette loop, which appears at position 35. 



Order 


Specification 




(72[l + ...+i] 


1 


{1123431142} 


8.7 (±4.4) 


8.7 (±4.4) 


2 


{111222111222} 


82.7 (±9.2) 


7.2 (±4.1) 


3 


{123413143213} 


12.6 (±4.4) 


6.0 (±4.0) 


4 


{123123} 


14.4 (±5.0) 


5.8 (±3.9) 


5 


{12131213} 


9.3 (±4.5) 


5.6 (±3.9) 


6 


{12123434} 


22.5 (±5.6) 


5.3 (±3.8) 


7 


{11121112} 


19.0 (±4.6) 


5.1 (±3.7) 


8 


{112112334334} 


17.1 (±4.8) 


5.0 (±3.8) 


9 


{12332133} 


9.5 (±4.6) 


4.8 (±3.7) 


10 


{111122111122} 


60.1 (±8.4) 


4.6 (±3.6) 


11 


{112334331142} 


9.6 (±4.6) 


4.4 (±3.7) 


12 


{1122322113} 


14.4 (±4.6) 


4.0 (±3.9) 


13 


{112112} 


15.2 (±4.5) 


3.6 (±3.7) 


14 


{12121313} 


50.1 (±7.5) 


3.5 (±3.8) 


15 


{12131213} 


20.0 (±5.2) 


3.2 (±3.9) 


16 


{1221221313} 


33.1 (±5.6) 


3.1 (±3.8) 


17 


{1123311233} 


17.7 (±4.6) 


3.0 (±3.8) 


18 


{112233221133} 


23.0 (±4.7) 


2.7 (±3.8) 


19 


{1123321133} 


14.9 (±4.6) 


2.6 (±3.7) 


20 


{122133122133} 


24.3 (±5.2) 


2.6 (±3.7) 


21 


{123213} 


14.2 (±5.1) 


2.5 (±3.6) 


22 


{1213312133} 


14.5 (±4.4) 


2.5 (±3.7) 


23 


{12343142} 


12.6 (±4.8) 


2.4 (±3.8) 


24 


{12121212} 


52.8 (±7.8) 


2.4 (±3.8) 


25 


{12343214} 


12.6 (±4.3) 


2.4 (±3.8) 


26 


{1221312213} 


10.7 (±4.4) 


2.3 (±3.7) 


27 


{1112211122} 


32.6 (±5.6) 


2.3 (±3.8) 


28 


{11231123} 


10.5 (±4.7) 


2.3 (±3.7) 


29 


{112334332114} 


12.3 (±4.5) 


2.3 (±3.8) 


30 


{1221312213} 


10.9 (±4.5) 


2.3 (±3.8) 


31 


{11232113} 


10.9 (±4.4) 


2.3 (±4.0) 


32 


{1123432114} 


10.4 (±4.4) 


2.3 (±4.0) 


33 


{11221122} 


26.6 (±4.9) 


2.2 (±4.0) 


34 


{12121212} 


71.4 (±9.4) 


2.2 (±3.9) 


35 


{1212} 


22.6 (±5.6) 


2.2 (±4.0) 




T 


295.0 (±20.0) 
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FIG. 1. Effective meson masses in the Schwinger model calculated with a true action (full 
dynamical fermion) simulation (/3 = 2.29, k = 0.26) and with a quenched approximate action 
simulation at /? = 2.50. 
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• Full fermion 




FIG. 2. The static potential in physical units calculated in three different ways. Crosses 
indicate the quenched, /? = 2.50 result, open circles indicate the deficit corrected quenched results 
and the filled circles are the full dynamical fermion results for /3 = 2.29, k = 0.25. The curves are 
described in the text. 
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Loop 1 



Loop 2 



Loop 3 



Loop 3' 









FIG. 3. Loops included in the first attempt to improve on the quenched approximation for the 
Schwinger model. 
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FIG. 4. The static potential in lattice units using different orders of approximation. Open 
symbols denote raw measurements at a given order, while full symbols denote deficit corrected 
results. The corresponding full dynamical simulation results are indicated by the crosses at the 
right of the diagram. 
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1st: 3x3 



3rd: 5x5 





2nd: 1 x 6 



4th: 1 x2 



FIG. 5. Loops which give the optimum approximation to the true Schwinger model dynamical 
fermion action. 
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